Derivation on Optical bench matrix system.

Photon polariazation jones matrices.



In [1]:

    
%pylab inline
from sympy import MatrixSymbol, Matrix, Symbol
from sympy import S, simplify, count_ops, oo
from sympy.physics.quantum import TensorProduct
from sympy import *

import sys
sys.path.append('../')
import Icarus

or2 = 1.0/sqrt(2)









    



Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.



In [2]:

    
def mapping_matrix_test(cases, U):
    """
        Tests that for each case the map works.
    """
    
    for case in cases:
        a = case[0]
        b = case[1]
        
        assert U*a == b, "Test failed"

def mapping_matrix(cases):
    """
        Create matrix U such that U*A = B.
    """
    
    l = len(cases[0][0])

    for case in cases:
        assert len(case[0]) == l
        assert len(case[1]) == l        
    
    for k in xrange(len(cases)):
    
        a = cases[k][0]
        b = cases[k][1]
        
        if k == 0:
            U = b * a.T
        else:
            U += b * a.T

    mapping_matrix_test(cases, U)
    return U



In [3]:

    
# Polarizations.

H = Matrix([[1], [0]])
V = Matrix([[0], [1]])
D = or2*(H+V)
A = or2*(H-V)
R = or2*(H+1j*V)
L = or2*(H-1j*V)

# Directions

i = Matrix([[1], [0]])
j = Matrix([[0], [1]])

# Energy

X = Matrix([[1], [0]])
XX = Matrix([[0], [1]])

# Null

null_matrix = Matrix([[0],[0]])


# Extended kronecker delta states

ih = TensorProduct(i, H)
iv = TensorProduct(i, V)
jh = TensorProduct(j, H)
jv = TensorProduct(j, V)

ix = TensorProduct(i, X)
ixx = TensorProduct(i, XX)
jx = TensorProduct(j, X)
jxx = TensorProduct(j, XX)


ixh = TensorProduct(X, ih)
ixxh = TensorProduct(XX, ih)
ixv = TensorProduct(X, iv)
ixxv = TensorProduct(XX, iv)

jxh = TensorProduct(X, jh)
jxxh = TensorProduct(XX, jh)
jxv = TensorProduct(X, jv)
jxxv = TensorProduct(XX, jv)

null_state = TensorProduct(null_matrix, TensorProduct(null_matrix, null_matrix))

Half and quarter wave plates.



In [4]:

    
def h(a = 0.0):
    return Matrix([
        [cos(2*a), sin(2*a)],
        [sin(2*a), - cos(2*a)]
    ])

def q(a = 0.0):
    return or2*Matrix([
        [1.0 + I*cos(2*a), I*sin(2*a)],
        [1.0j*sin(2*a), 1.0 - I*cos(2*a)],
    ])



In [5]:

    
def test_hwp():
    assert h(pi/8)*H - D == Matrix([[0],[0]])
    assert h(pi/8)*V - A == Matrix([[0],[0]])
    
def test_qwp():
    assert q(pi/4)*H - R == Matrix([[0],[0]])
    assert q(pi/4)*V - 1j*L == Matrix([[0],[0]])



In [6]:

    
test_hwp()
test_qwp()

alpha1 = Symbol('alpha1')
alpha2 = Symbol('alpha2')
beta1 = Symbol('beta1')
beta2 = Symbol('beta2')

ah = TensorProduct(X, TensorProduct(i, h(alpha1))) 
bh = TensorProduct(X, TensorProduct(j, h(alpha1)))
ch = TensorProduct(XX, TensorProduct(i, h(alpha2)))
dh = TensorProduct(XX, TensorProduct(j, h(alpha2)))

aq = TensorProduct(X, TensorProduct(i, q(beta1))) 
bq = TensorProduct(X, TensorProduct(j, q(beta1)))
cq = TensorProduct(XX, TensorProduct(i, q(beta2)))
dq = TensorProduct(XX, TensorProduct(j, q(beta2)))



In [7]:

    
zh = Matrix(8,8, np.zeros(64))
zh[:,0:2] += ah
zh[:,2:4] += bh
zh[:,4:6] += ch
zh[:,6:8] += dh


zq = Matrix(8,8, np.zeros(64))
zq[:,0:2] += aq
zq[:,2:4] += bq
zq[:,4:6] += cq
zq[:,6:8] += dq


HWP = zh
QWP = zq

Neutral beam splitter.

Does nothing but send half the photon state in the $\hat{i}$ direction and half in the $\hat{j}$ direction.

$ \langle NBS \ | \ i \ \psi \rangle = \frac{1}{\sqrt{2}}( | \ i \ \psi \rangle + | \ j \ \psi \rangle )$

and

$ \langle NBS \ | \ j \ \psi \rangle = \frac{1}{\sqrt{2}}( | \ i \ \psi \rangle + | \ j \ \psi \rangle )$



In [8]:

    
NBS_cases = [
    [ixh, or2*(ixh + jxh)],
    [jxh, or2*(ixh + jxh)],
    [ixxh, or2*(ixxh + jxxh)],
    [jxxh, or2*(ixxh + jxxh)],
    
    [ixv, or2*(ixv + jxv)],
    [jxv, or2*(ixv + jxv)],
    [ixxv, or2*(ixxv + jxxv)],
    [jxxv, or2*(ixxv + jxxv)],
]

NBS = mapping_matrix(NBS_cases)

Spectrometers.

One channel allows excitons to pass, but blocks biexcitons. The other channel allows biexcitons to pass and blocks excitons.



In [9]:

    
Si_cases = [
    [ixxh, ixxh],
    [ixxv, ixxv],
    [ixh, null_state],
    [ixv, null_state],
]

Si = mapping_matrix(Si_cases)

Sj_cases = [
    [jxh, jxh],
    [jxv, jxv],
    [jxxh, null_state],
    [jxxv, null_state],
]

Si = mapping_matrix(Si_cases)
Sj = mapping_matrix(Sj_cases)

S = Si + Sj

Polarising beam splitter.

Changes to photon path depending on its polarization.



In [10]:

    
PBSi_cases = [
    [ixxh, ixxh],
    [ixxv, jxxv],
    
    [ixh, ixh],
    [ixv, jxv],
]

PBSj_cases = [
    [jxh, jxh],
    [jxv, ixv],
    [jxxh, jxxh],
    [jxxv, ixxv],

]

PBSi = mapping_matrix(PBSi_cases)
PBSj = mapping_matrix(PBSj_cases)

PBS = PBSi + PBSj

Biphoton states



In [11]:

    
HWPHWP = TensorProduct(HWP, HWP)
QWPQWP = TensorProduct(QWP, QWP)
NBSNBS = TensorProduct(NBS, NBS)
SS = TensorProduct(S, S)
PBSPBS = TensorProduct(PBS, PBS)



In [12]:

    
phi = Symbol('phi')
state = or2*(TensorProduct(ixh, ixxh) + phi* TensorProduct(ixv, ixxv))



In [13]:

    
D1 = jxh
D2 = ixv
D3 = ixxh
D4 = jxxv

D1D3 = TensorProduct(D1, D3)
D1D4 = TensorProduct(D1, D4)
D2D3 = TensorProduct(D2, D3)
D2D4 = TensorProduct(D2, D4)



In [14]:

    
p_state = PBSPBS*SS*QWPQWP*HWPHWP*NBSNBS*state



In [15]:

    
d1d3p = (4*abs((D1D3.T * p_state)[0])**2)
d1d4p = (4*abs((D1D4.T * p_state)[0])**2)
d2d3p = (4*abs((D2D3.T * p_state)[0])**2)
d2d4p = (4*abs((D2D4.T * p_state)[0])**2)



In [16]:

    
d1d3e = lambdify((phi, alpha1, alpha2, beta1, beta2), d1d3p, 'numpy')
d1d4e = lambdify((phi, alpha1, alpha2, beta1, beta2), d1d4p, 'numpy')
d2d3e = lambdify((phi, alpha1, alpha2, beta1, beta2), d2d3p, 'numpy')
d2d4e = lambdify((phi, alpha1, alpha2, beta1, beta2), d2d4p, 'numpy')



In [17]:

    
def foo(phase, hwpx, hwpxx, qwpx, qwpxx):
    
    p1 = d1d3e(phase, hwpx, hwpxx, qwpx, qwpxx)
    p2 = d1d4e(phase, hwpx, hwpxx, qwpx, qwpxx)
    p3 = d2d3e(phase, hwpx, hwpxx, qwpx, qwpxx)
    p4 = d2d4e(phase, hwpx, hwpxx, qwpx, qwpxx)
    
    return np.array([p1, p2, p3, p4])

def bar(phase, hwpx, hwpxx, qwpx, qwpxx):
    
    a1, a2, a3, a4 = foo(phase, hwpx, hwpxx, qwpx, qwpxx)
    
    return (a1-a2)/(a1+a2), (a3-a4)/(a3+a4)



In [19]:

    
measurements = [
   [0, 0, 0, 0],
   [0, 0, np.pi/4, 0],
   [np.pi/4, 0, np.pi/4, 0],
   [np.pi/4, 0, 0, 0],
   [0, np.pi/4, 0, 0],
   [0, np.pi/4, np.pi/4, 0],
   [np.pi/8, np.pi/4, np.pi/4, 0],
   [np.pi/8, np.pi/4, 0, 0],
   [np.pi/8, np.pi/4, 0, np.pi/4],
   [np.pi/8, np.pi/4, np.pi/8, np.pi/4],
   [0, np.pi/4, np.pi/8, np.pi/4],
   [0, 0, np.pi/8, np.pi/4],
   [np.pi/4, 0, np.pi/8, np.pi/4],
   [np.pi/4, 0, 0, -np.pi/4],
   [0, 0, 0, -np.pi/4],
   [0, np.pi/4, 0, -np.pi/4]
   ]



In [20]:

    
angles = np.linspace(-np.pi, np.pi, 100)

plt.plot([bar(1, m[0], m[2], m[1], m[3]) for m in measurements])

plt.ylim([-1.05, 1.05])
plt.legend(['deg1', 'deg2'])









    Out[20]:





<matplotlib.legend.Legend at 0x10707be10>



In [114]:

    
t = np.linspace(0, 10, 100)

qd = Icarus.QuantumDot()

def make_phase(fss, crosstau = 1e10):
    qd.FSS = fss
    qd.xlifetime = 1.
    qd.crosstau = crosstau
    return qd.generate_phase()

phases = np.array([make_phase(0, 3) for p in xrange(100)])



In [115]:

    
angles = np.linspace(-np.pi, np.pi, 200)

plt.plot(angles, [bar(make_phase(0, 1e10), np.pi/8, np.pi/4, np.pi, angle)[0] for angle in angles])
plt.plot(angles, [bar(phases.mean(), np.pi/8, np.pi/4, np.pi, angle)[0] for angle in angles])

plt.ylim([-1.05, 1.05])
plt.legend(['no cross', 'with_cross'])









    Out[115]:





<matplotlib.legend.Legend at 0x10ad05d50>



In [ ]: